Bistability signatures in nonequilibrium charge transport through molecular quantum 

dots 
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We investigate the transient nonequilibrium dynamics of a molecular junction biased by a fi- 
nite voltage and strongly coupled to internal vibrational degrees of freedom. Using two different, 
numerical exact techniques, diagrammatic Monte Carlo and the multilayer multiconfiguration time- 
dependent Hartree method, we show that the steady state current through the junction may depend 
sensitively on the initial preparation of the system, thus revealing signatures of bistability. The in- 
fluence of the bias voltage and the transient dynamics on the phenomenon of bistability is analyzed. 
Furthermore, a possible relation to the phenomenon of stochastic switching in nanoelectromecanical 
devices is discussed. 

PACS numbers: 85.65.+h, 73.63.-b, 63.22.-m 



Since the first realization of a single-molecule junction 
the field of molecular electronics has seen a rapid devel- 
opment [l|, Q. Experimental investigations of the con- 
ductance properties of single-molecule junctions have re- 
vealed a wealth of intriguing transport phenomena 0, 0] . 

Molecules can be considered as very small quantum 
dots. An important aspect that distinguishes them from 
semiconductor-based mesoscopic systems is the influence 
of the vibrational degrees of freedom. Due to the small 
size of molecules, the charging of the molecular bridge 
is often accompanied by a significant change of the nu- 
clear geometry, which indicates strong coupling between 
electronic and vibrational degrees of freedom. It mani- 
fests itself in the current-voltage characteristics of molec- 
ular junctions 0,(1] @ and may result in current-induced 
heating of the molecular bridge Q and large fluctuations 
Conformational changes of the geometry of the con- 
ducting molecule are possible mechanisms for switching 
behavior and negative differential resistance 

The experimental findings have stimulated great inter- 
est in the basic mechanisms of quantum transport at the 
molecular scale, in particular effects due to electronic- 
vibrational coupling. A variety of different theoretical 
methods have been applied to study these phenomena, 
including scattering theory, nonequilibrium Green's func- 
tion approaches, and master equation methods (see, e. g., 
and references therein) . Although much physical in- 
sight has been obtained by the application of these meth- 
ods, all these approaches involve significant approxima- 
tions. 

To elucidate the detailed mechanisms and address the 
full complexity of the nonequilibrium transport prob- 
lem, advanced numerical techniques that do not in- 
volve intrinsic approximation are required. Powerful 



methods that have been proposed in this context are 
the diagrammatic Monte Carlo simulation (diagMC) ap- 
proach [9( , multilayer multiconfiguration time-dependent 
Hartree method in second quantization representation 
(ML-MCTDH-SQR) 0, [Tl| as well as standard quan- 



tum Monte Carlo and iterative path integral summation 
scheme 12, HI- 

Here, we employ the first two methods to address 
two important and largely unsolved questions: (i) How 
is the steady state in a molecular quantum dot estab- 
lished starting from a specified preparation? (ii) Does 
the steady state depend on the initial preparation, e.g. 
on the initial occupation of the dot? In particular (ii), 
which is closely related to the phenomena of bistability 
and hysteresis, has been discussed controversially based 
on approximate methods 3"l3|- On the other hand, 
there are a number of works in which hysteretic behavior 
was observed experimentally Our findings indicate 
that, for certain parameter regimes, a remarkable mem- 
ory effect in the nonequilibrium dynamics exists, which 
manifests itself in different steady states for different ini- 
tial preparations. We discuss the relation of this finding 
to earlier predictions of bistability for the model inves- 
tigated in 14- HI as well as to stochastic switching ob- 
served in nanoelectromechanical systems fl9j - [23| . 



In order to study vibrationally coupled electron trans- 
port in a molecular quantum dot we employ the resonant 
tunneling model given by the Hamiltonian [H, [13, HH , 
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which describes a single electronic level of energy cd 
(we shall consider a spinless system and use units in 
which h = e = ks = I throughout), Hr> = eod^d, cou- 
pled to two noninteracting fermionic reservoirs, Hlr = 
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J2 a .k e aka i ak a ak , describing the left and right (a = L,R) 
metallic electrodes with energy dispersion e a k, and a 
bosonic part, i? p h = XL ^^bi which models the vibra- 
tional degrees of freedom of the molecule within the har- 
monic approximation employing normal modes with fre- 
quencies lj k . The fermionic environments serve as charge 
reservoirs, inducing a nonequilibrium current by virtue of 
the coupling between the leads and the molecular quan- 
tum dot 
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where t a k denotes the dot-lead coupling strength between 
the fcth electronic level of lead a and the molecule. Fi- 
nally, electronic-vibrational coupling is described by the 
interaction part 



(3) 



with coupling constants X K . The last term is a 
counter term that corresponds to the polaron shift 
^ K Aj/w K . The molecule- lead and electron-vibrational 
couplings are characterized by the hybridization L Q (e) = 
2 71 " J2k \^ak\ 2 S(e — e a k) and spectral density J{ui) = 
7r J2 K ^k^( w ~ w «)> respectively. In our numerical studies 
we use different models for the molecule-lead coupling; 
on the qualitative level the phenomena we observe are 
universal and independent of them. 

To characterize the nonequilibrium dynamics of the 
quantum dot, we consider the dot population, P(t) = 
{(fi(t)d{t)) = ti{po<tf(t)d(t)}, and the electrical current 
through the dot, I{t) = (l/2){d/dt)(N L (t) - N R (t)) = 
[I L (t) - I R (t)]/2. Thereby, N a = Y,k a tk a ak measures 
the number of particles in contact a, I a (t) are the cur- 
rents from the individual electrodes to the dot, and the 
initial preparation is given by the density operator po 
at time t = 0, which describes a factorizing initial state, 
with the dot being either empty or occupied. Initially the 
leads are in thermal equilibrium with energies shifted ac- 
cording to the chemical potentials fj, a ; alternatively, one 
could keep each of the leads filled up to different chem- 
ical potentials. In the wide band limit, the stationary 
currents from these two initial conditions have negligible 
numerical differences (2||. The bias V = /i£ — /j,r en- 
sures that for sufficiently long times a stationary state is 
reached not only for the average current I(t) , but also 
for the currents from the individual electrodes to the dot 
I a (t) , which are equal in the limit t —¥ oo [3l| . In practice 
I(t) converges to the steady state value much faster than 
I a (t) though, see Fig. Q] and [l,[loj]. 

The reached stationary nonequilibrium state is, if er- 
godicity holds, supposed to be unique and independent 
of the initial preparation of the system. Therefore, even 
though different preparations might yield different tran- 
sient dynamics, in the long-time limit one expects to find 



the transport properties independent of the evolution his- 
tory. In this Letter we investigate the validity of this 
conjecture. 

We employ two different numerical approaches, the 
ML-MCTDH-SQR method and the diagMC approach 
(see for details). Briefly, the ML-MCTDH-SQR 

method is a variational basis-set approach to study quan- 
tum dynamics for large systems containing identical par- 
ticles. Within this approach the wave function is repre- 
sented by a recursive, layered expansion in Fock space 
employing the occupation number basis. Its time evo- 
lution is then determined from the Dirac-Frenkel vari- 
ational principle by dynamically optimizing all the pa- 
rameters. On the other hand, the diagMC method relies 
on an expansion of the time evolution in terms of the 
dot-lead coupling. After integrating out all environmen- 
tal degrees of freedom, one obtains an infinite sum over 
Feynman diagrams which is then evaluated by a stochas- 
tic MC scheme. Both approaches are numerically exact, 
as the respective errors can be made arbitrarily small: 
for diagMC, the statistical error due to the MC sam- 
pling can be easily minimized by increasing the number 
of MC measurements, while ML-MCTDH-SQR results 
can be systematically converged by increasing the num- 
ber of states as well as the basis set. Both approaches 
yield consistent results for the transient dynamics as well 
as for the stationary state. Fig. [T] depicts results of both 
approaches for the time-dependent current for a model 
with a semielliptic molecule-lead coupling and an Ohmic 
spectral density, showing excellent agreement. 

While the results of Fig. [1] corresponding to off- 
resonant transport, are an example for a stationary state 
that is independent of the initial preparation, various 
claims have been made regarding the existence of bista- 
bility or hysteresis in the model under investigation [hj - 



llq . |27I |. These claims, however, received considerable 
criticism; the notion of bistable regimes in this model 
is anything but generally accepted |28| ■ Usually, the rea- 
soning in favor of such effects assumes the existence of 
a parameter regime (typically requiring strong electron- 
phonon coupling) within which the effective potential of 
the phonon mode becomes multistable even in the sta- 
tionary regime. Accordingly, upon entering this regime 
from the smaller or larger voltage domain, the system 
ends up in different potential minima, resulting in differ- 
ent stationary states and corresponding different trans- 
port properties for otherwise identical parameters 14j . 

These claims have been debated for almost a decade 
without reaching an accepted conclusion. On the one 
side, a general proof of non-existence of bistability is 
hard to give, while on the other side all previous works 
in favor of the phenomenon rely on approximative meth- 
ods (e. g. the Born-Oppenheimer approximation of [Hj]) 
whose accuracy is difficult to judge. With our numeri- 
cally exact approaches at hand, however, we can give a 
more detailed account of the underlying dynamics in the 
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FIG. 1: (color online) Comparison of ML-MCTDH-SQR 
(lines) and diagMC data (symbols) for an initially empty 
(black) and occupied (red/grey) dot for I(t) (circles/solid 
lines), (squares/dashed lines), and —In(t) (dia- 

monds/dotted lines). The main graph shows the average 
current for e D = 4.7T, V = 1.25r, F a (e) = - (e/e c ) 2 

for |e| < e c with bandwidth 2e c = 50r (semielliptic band), 
and an Ohmic spectral density J(cj) = 2nu) exp(— uj/uj c ) with 
u)c = 0.775F at zero temperature; filled and open symbols 
refer to I(t) and [Ihit) — In(t)]/2, respectively. The inset dis- 
plays the different timescales on which I(t), Ihit), and In(t) 
reach the stationary regime, while the upper panel demon- 
strates that all currents converge to the same stationary value. 



proposed regimes of bistability. According to [14J the two 
points in parameter space which minimize the action and 
result in two different values for the current correspond 
to a nearly empty and a nearly fully occupied dot. Thus, 
numerical simulations with initial occupation numbers 
and 1 should lead to the different dot configurations re- 
quired for bistability to occur. Therefore, in the following 
we investigate the influence of the initial preparation on 
the stationary state, which, in the absence of bistability, 
should not exist. Since, as noted above, in the longtime 
limit all currents I a (t), I(t) converge to the same value, 
here we restrict ourselves to the stationary behavior of 
/(*). 

We consider a model with a single, strongly-coupled 
vibrational mode, J(oj) = itXq6(oj — uj ). Fig. [2] shows 
the time-dependent current for the two different initial 
preparations (i.e. empty and occupied dot) at different 
voltages for parameters in the nonadiabatic regime, i. e. 
cjo > r. Strikingly, the currents converge towards sta- 
tionary values on timescales roughly of the order of 
for almost all voltages. This is in accordance with many 
previous studies of the real-time dynamics, both for the 
as well as for the Anderson model 



present [29( as well as for the Anderson model [26[ and 
can be rationalized by noting that the strength of the tun- 
neling coupling T defines the typical timescale for charge 
dynamics. However, except for the high voltages the sta- 
tionary values of the current are clearly different, with 
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FIG. 2: Current I{t) for e D = -ST, T a (e) = V/2 for |e| < e c 
(fiat band), and a single phonon mode with Ao = 8r and 
cjo = 4r, for an initially empty (black) and occupied (red) 
dot at zero temperature. Open symbols refer to diagMC data 
for a smooth switch-on of the tunneling coupling, T(e,t) = 
sin 2 [7rf/(2T sw )]r(e) for t < r sw (= Sr^ 1 for V = 5F and 
35r, else 6F _1 ), while filled symbols refer to an instantaneous 
switch-on. The bandwidth is 2e c = 41 F for V = 35F, else 26F. 



the initially empty dot leading to a significantly larger 
current than the initially occupied one. This intriguing 
finding contradicts the notion of an unique stationary 
state and instead supports the idea of bistability. A sim- 
ilar result is obtained in the adiabatic regime of rather 
slow vibrational motion, i.e. ujq < T, depicted in Fig|3] 
This is the regime originally suggested for the existence 
of bistability in the theoretical approaches 14| . The nu- 
merical results in Figs JUGS indicate the dependence of the 
long-time current on the initial occupation of the dot. 
This finding is corroborated by the dynamics of the pop- 
ulation of the dot (data not shown) . The population of 
the dot may converge to its stationary value on a signif- 
icantly longer timescale than the current. However, for 
the examples considered above its stationary state ex- 
hibits a dependence on the initial occupation even more 
pronounced than for the current. 

The comparison of the different panels in Fig. [5] shows 
that the sensitivity of the dependence of the current on 
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FIG. 3: Similar to Fig. d but with ML-MCTDH-SQR result 
for the parameters: e_o = — 2.5r, uia = T/8, Ao = T, and a 
bandwidth 2e c = 20r. Insets: Magnified versions of the same 
plots in the stationary regime. 



the initial preparation is influenced by the applied bias 
voltage. While for V = 35T both currents are approach- 
ing each other within a timescale r st of the order of T _1 
they still seem to be transient for V = 20r within the 
times accessible by our simulations. Here, it is not clear 
whether both currents converge to different or the same 
stationary value. In the latter case the corresponding 
stationary timescale r st would be much larger than T -1 , 
offering an alternative explanation for the mismatch in 
the current for not too large voltages in both, Figs. [2] 
and [3] Instead of identifying the plateau value as the 
actual steady state, it could be interpreted as a transient 
current with a timescale r s t exceeding T -1 by orders of 
magnitude. 

The temperature T of the leads has a similar influence 
on the effect: an increase of T destroys the bistability sig- 
natures just as an increase of V does. On the other hand, 
it is noted that the difference in the long-time behavior is 
independent of the details of the initial switch-on of the 
tunneling coupling (cf. Fig. [2]) as well as of the specific 
form of the hybridization (i.e. semielliptic vs. flat band). 

We would like to stress that the numerically exact, 
time-dependent methods employed in our studies cannot 
directly address the stationary state, which is formally 



obtained in the infinite time limit, but only map the dy- 
namics over long but finite timescales. It is thus not pos- 
sible to unambiguously determine whether the differences 
in the long-time transport behavior account for truly 
distinct stationary states, or whether the corresponding 
transport properties would still decay to the same unique 
stationary value, void of any imprints of the preparation, 
on timescales significantly longer than those covered by 
our calculations. While in the latter case the concept of 
bistability does not have to be introduced, an explana- 
tion for this surprisingly strong separation of timescales 
is, to our knowledge, still lacking - that is, the timescales 
necessary to reach the true stationary state within and 
outside the aforementioned parameter regime. 

Since the pronounced transient dynamics always die 
out on the same timescale (of the order of T -1 ), we con- 
clude that the effect is not caused but merely triggered by 
the changes in the initial preparations. Furthermore, the 
phenomenon seems to be completely unaffected by details 
of the transient dynamics. For voltages at which we ob- 
serve the separation of timescales, the current originating 
from an instantaneous initial switch-on of the tunneling 
coupling exhibits very rich transient dynamics, suggest- 
ing a rather extensive exploration of the effective poten- 
tial surfaces, while for a continuous switch-on process, 
one finds quite a smooth convergence process. Yet in 
both cases we observe the same separation of timescales. 
Since the system lacks any energy scale which can be set 
in relation to this novel long timescale it should, if finite, 
be generated in a yet unknown way. 

We also would like to mention that the described phe- 
nomenon is not directly related to the stochastic switch- 
ing processes leading among other things to random tele- 
graph noise frequently discussed in the context of shut- 



tling in nanoelectromechanical devices 



While 

the switching times would be of the order of 10 3 — 10 4 in 
units of r _1 for most of the plots shown above (we use 
the procedure from [2l|) and thus apparently in accor- 
dance with the reported data, we would like to point out 
that our simulations would then describe the very onset 
of the stochastic switching process, which is not covered 
by any of the existing studies. 

We finally comment on the applicability of our findings, 
obtained for the standard model of vibrationally-coupled 
electron transport, to realistic molecular junctions. Typ- 
ical values of the coupling T vary between meV and a few 
eV, depending on the specifics of the binding group and 
the geometry. Assuming an average value of T = 0.1 eV 
(which is realized, e. g. in benzenealkenethiol-gold junc- 
tions [131), our results predict the phenomenon of bista- 
bility to occur both in the nonadiabatic regime (e. g. for 
a vibrational mode with a higher frequency of ujq = 0.4 
eV as in Fig. [2]) as well as in the adiabatic regime (e. g. 
uiq = 0.0125 eV as in Fig. [3]) for sufficiently large elec- 
tronic vibrational coupling {X/ujq > 1) and not too large 
voltages (V < 2 V). These parameters are expected to be 
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of relevance for many molecular junctions, e.g. for some 
of the molecular j unctions ,where hysteretic behavior was 
observed experimentally [181 ] . According to our results, 
the phenomenon of bistability persists on a timescale of 
at least a picosecond, which is longer than typical vibra- 
tional periods. It should be observable with methods of 
femtosecond spectroscopy. In view of the current exper- 
imental progress (zj, spectrosopic techniques that allow 
to directly monitor ultrafast processes in molecular junc- 
tions are expected to become available in the near future. 

To summarize, we have presented numerically exact 
results for the nonequilibrium dynamics of a molecular 
quantum dot. In a wide parameter range we find a dis- 
tinct dependence of the steady state current on the ini- 
tial preparation of the molecular junction. Our analysis 
shows that this phenomenon is related to the previously 
predicted bistability of the model. In contrast to the 
previous approximate approaches, the present numeri- 
cally exact results unambiguously prove the existence of 
bistability signatures in the dynamics over a time scale 
significantly longer than the tunneling time. 
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